#!/usr/bin/env perl

#test the mapping with two hexagonal ring structures, refer to 'hexagonalMapping.txt' in 'testData'
#for a listing of all valid mappings. In a second step, an additional domain is added, such that
#symmetry is broken and only one unique mapping exists.

use QueryProcessor;
use SubcomplexMerger;
use SubcomplexFinder;
use Scorer;

use SBG::Domain::Sphere;
use SBG::DomainIO::cofm;

use Test::More;
use Set::Scalar;
use PDL;
use PDL::Transform;
use File::Basename qw/dirname/;
use File::Spec::Functions;
use constant PI => 3.14159265;

use strict;
use warnings;

#'inf' not accepted by 'strict'
my $inf = 1000000000;

#radius, chain ids
my $r = 20;
my @chains=('A','B','C','D','E','F');

#z-offset, rotation angles
my $z1 = 20;
my $z2 = -20;
my $phi1x = 20;
my $phi1y = -7;
my $phi1z = 40;
my $phi2x = 10;
my $phi2y = 38;
my $phi2z = 45;

my $nDomains = scalar(@chains);

#create first hexagonal ring + 1 domain that is added later
my @doms1 = ();
for(my $i=0; $i<$nDomains; $i++){
	my $x = $r*cos($i*2*PI/$nDomains);
	my $y = $r*sin($i*2*PI/$nDomains);
	my $dom = SBG::Domain::Sphere->new();
	my $c1= PDL->new([$x,$y,$z1,1]);
	$c1 = $c1->transpose();
	
	#apply transformation
	my $Tx = pdl([[1,0,0,0],[0,cos($phi1x),-sin($phi1x),0],[0,sin($phi1x),cos($phi1x),0],[0,0,0,1]]);
	my $Ty = pdl([[cos($phi1y),0,sin($phi1y),0],[0,1,0,0],[-sin($phi1y),0,cos($phi1y),0],[0,0,0,1]]);
	my $Tz = pdl([[cos($phi1z),-sin($phi1z),0,0],[sin($phi1z),cos($phi1z),0,0],[0,0,1,0],[0,0,0,1]]);
	
	my $c2 = $Tz x $Ty x $Tx x $c1;
	$c2 = $c2->transpose();
	
	#as crosshairs are not used here: take the centroid coordinates for all of the 7 atoms
	my $coords = $c2->glue(1, $c2, $c2, $c2, $c2, $c2, $c2);
	$dom->coords($coords);
	my $chain = $chains[$i];
	$dom->descriptor("CHAIN $chain");
	$dom->classification('a');
	$dom->model(1);
	$dom->assembly(1);
	$dom->pdbid("1tes");	
	push @doms1, $dom;
}
#create additional domain 
my $additionalDom1 = SBG::Domain::Sphere->new();
my $d1 = $doms1[1]->coords - $doms1[0]->coords;
$additionalDom1->coords($doms1[0]->coords+$d1*2/3);
$additionalDom1->descriptor("CHAIN G");
$additionalDom1->classification('a');
$additionalDom1->model(1);
$additionalDom1->assembly(1);
$additionalDom1->pdbid("1tes");	

#create second hexagonal ring + 1 domain that is added later
my @doms2 = ();
my $additional2 = SBG::Domain::Sphere->new();
for(my $i=0; $i<$nDomains; $i++){
	my $x = $r*cos($i*2*PI/$nDomains);
	my $y = $r*sin($i*2*PI/$nDomains);
	my $dom = SBG::Domain::Sphere->new();
	my $c1= PDL->new([$x,$y,$z1,1]);
	$c1 = $c1->transpose();

	#apply transformation
	my $Tx = pdl([[1,0,0,0],[0,cos($phi2x),-sin($phi2x),0],[0,sin($phi2x),cos($phi2x),0],[0,0,0,1]]);
	my $Ty = pdl([[cos($phi2y),0,sin($phi2y),0],[0,1,0,0],[-sin($phi2y),0,cos($phi2y),0],[0,0,0,1]]);
	my $Tz = pdl([[cos($phi2z),-sin($phi2z),0,0],[sin($phi2z),cos($phi2z),0,0],[0,0,1,0],[0,0,0,1]]);
	
	my $c2 = $Tz x $Ty x $Tx x $c1;
	$c2 = $c2->transpose();
	
	my $crosshair = PDL->new([0,0,0,1],[0,0,0,1],[0,0,0,1],[0,0,0,1],[0,0,0,1],[0,0,0,1]);
	my $coords = $c2->glue(1, $c2, $c2, $c2, $c2, $c2, $c2);
	$dom->coords($coords);
	my $chain = $chains[$i];
	$dom->descriptor("CHAIN $chain");
	$dom->classification('a');
	$dom->model(1);
	$dom->assembly(1);
	$dom->pdbid("2tes");
	push @doms2, $dom;
}
#create additional domain 
my $additionalDom2 = SBG::Domain::Sphere->new();
my $d2 = $doms2[1]->coords - $doms2[0]->coords;
$additionalDom2->coords($doms2[0]->coords+$d2*2/3);
$additionalDom2->descriptor("CHAIN G");
$additionalDom2->classification('a');
$additionalDom2->model(1);
$additionalDom2->assembly(1);
$additionalDom2->pdbid("1tes");	

#initilialize modules
my $scorer = Scorer->new(maxRMSD=>5);
my $subcomplexMerger = SubcomplexMerger->new(scorer=>$scorer);
my $subcomplexFinder = SubcomplexFinder->new(scorer=>$scorer);

#get initial subcomplexes
my @subcomplexes = $subcomplexFinder->pairwise(\@doms1,\@doms2);
#merge till last stage
my $loop=1;
while($loop){
	my ($finished, @merged) = $subcomplexMerger->mergePrefixMethod(\@subcomplexes, $inf);
	if(scalar(@merged) != 0){
		@subcomplexes = @merged;
	}else{
		$loop=0;
	}
}

#extract detected similarities in string representation
my $detectedSubcomplexes=Set::Scalar->new();
for my $s (@subcomplexes){
	my $a='';
	my $b='';
	for my $domain (@{$s->domainsA}){
		$a= join '',$a, substr($domain->id,9,1);
	}
	for my $domain (@{$s->domainsB}){
		$b= join '',$b, substr($domain->id,9,1);
	}
	$detectedSubcomplexes->insert("$a-$b");
}

#check for false negatives and positives (only similarities with 6 components are considered)
my $correctSubcomplexes=Set::Scalar->new(
	"ABCDEF-ABCDEF",
	"ABCDEF-BCDEFA",
	"ABCDEF-CDEFAB",
	"ABCDEF-DEFABC",
	"ABCDEF-EFABCD",
	"ABCDEF-FABCDE",
	"ABCDEF-AFEDCB",
	"ABCDEF-BAFEDC",
	"ABCDEF-CBAFED",
	"ABCDEF-DCBAFE",
	"ABCDEF-EDCBAF",
	"ABCDEF-FEDCBA",
);

ok($detectedSubcomplexes == $correctSubcomplexes, "all correct subcomplexes identified");

#repeat the same, now with the the additional domain added, which should provide for a unique mapping
push @doms1, $additionalDom1;
push @doms2, $additionalDom2;

#get initial subcomplexes
my @subcomplexes2 = $subcomplexFinder->pairwise(\@doms1,\@doms2);
#merge till last stage
my $loop2=1;
while($loop2){
	my ($finished, @merged) = $subcomplexMerger->mergePrefixMethod(\@subcomplexes2, $inf);
	if(scalar(@merged) != 0){
		@subcomplexes2 = @merged;
	}else{
		$loop2=0;
	}
}

#extract detected similarities in string representation
my $detectedSubcomplexes2=Set::Scalar->new();
for my $s (@subcomplexes2){
	my $a='';
	my $b='';
	for my $domain (@{$s->domainsA}){
		$a= join '',$a, substr($domain->id,9,1);
	}
	for my $domain (@{$s->domainsB}){
		$b= join '',$b, substr($domain->id,9,1);
	}
	$detectedSubcomplexes2->insert("$a-$b");
}

#check for false negatives and positives (only similarities with 6 components are considered)
my $correctSubcomplexes2=Set::Scalar->new(
	"ABCDEFG-ABCDEFG"
);

ok($detectedSubcomplexes2 == $correctSubcomplexes2, "unique mapping correct");

done_testing;